      Subroutine RK4
      include "param.fi"
      include "common.fi"

      Integer I, ILevelRK4

      Double Precision kRho(NTOT3T, 4), oRho(NTOT3T)
      Double Precision kU1(NTOT3T, 4), oU1(NTOT3T)
      Double Precision kU2(NTOT3T, 4), oU2(NTOT3T)
      Double Precision kU3(NTOT3T, 4), oU3(NTOT3T)
      Double Precision kX1(NTOT3T, 4), oX1(NTOT3T)
      Double Precision kX2(NTOT3T, 4), oX2(NTOT3T)
      Double Precision kX3(NTOT3T, 4), oX3(NTOT3T)
      Double Precision dTempStep(4), dCoefStep(4)

      Do I = 1, NTOT3T
        oRho(I) = Rho(I)
        oU1(I) = U1(I)
        oU2(I) = U2(I)
        oU3(I) = U3(I)
        oX1(I) = X1(I)
        oX2(I) = X2(I)
        oX3(I) = X3(I)
      Enddo

      dTempStep(1) = 0.0
      dTempStep(2) = 0.5
      dTempStep(3) = 0.5
      dTempStep(4) = 1.0
      dCoefStep(1) = 1.0 / 6.0
      dCoefStep(2) = 1.0 / 3.0
      dCoefStep(3) = 1.0 / 3.0
      dCoefStep(4) = 1.0 / 6.0

      Do ILevelRK4 = 1, 4

      If (ILevelRK4 .GT. 1) Then
      Do I = 1, NTOT3T
        Rho(I) = oRho(I) 
     &         + dTempStep(ILevelRK4) * DT * kRho(I, ILevelRK4 - 1)
        If (Rho(I).LT.RDENSITYMIN) Rho(I) = RDENSITYMIN
        If (Rho(I).GT.RDENSITYMAX) Rho(I) = RDENSITYMAX
        U1(I) = oU1(I) 
     &         + dTempStep(ILevelRK4) * DT * kU1(I, ILevelRK4 - 1)
        U2(I) = oU2(I) 
     &         + dTempStep(ILevelRK4) * DT * kU2(I, ILevelRK4 - 1)
        U3(I) = oU3(I) 
     &         + dTempStep(ILevelRK4) * DT * kU3(I, ILevelRK4 - 1)
      Enddo

      If (ALE .EQ. ALELagrange) Then
        Do I = NTOT3T - NPWALL + 1, NTOT3T
          U1(I) = 0.0
          U2(I) = 0.0
          U3(I) = 0.0
        Enddo

      If (tagAVERMOVE .EQ. 0) Then
      Do I = 1, NTOT3T - NPWALL
        X1(I) = oX1(I)
     &         + dTempStep(ILevelRK4) * DT * kX1(I, ILevelRK4 - 1)
        If (X1(I).GT.XMAX) X1(I) = X1(I) - CUBESIZE
        If (X1(I).LT.XMIN) X1(I) = X1(I) + CUBESIZE
        If (DIMM.EQ.DIMM3D) Then
          X2(I) = oX2(I)
     &         + dTempStep(ILevelRK4) * DT * kX2(I, ILevelRK4 - 1)
          If (X2(I).GT.XMAX) X2(I) = X2(I) - CUBESIZE
          If (X2(I).LT.XMIN) X2(I) = X2(I) + CUBESIZE
        Endif
        X3(I) = oX3(I) 
     &         + dTempStep(ILevelRK4) * DT * kX3(I, ILevelRK4 - 1)
        If (X3(I).GT.XMAX) X3(I) = X3(I) - CUBESIZE
        If (X3(I).LT.XMIN) X3(I) = X3(I) + CUBESIZE
!        If (tagRegroup .NE. 0) Then
!          Call Regroup
!        Endif
      Enddo
      Else   ! tagAVERMOVE
      CALL AVER(kX1(1, ILevelRK4 - 1), TEMP1)
      CALL AVER(kX2(1, ILevelRK4 - 1), TEMP2)
      CALL AVER(kX3(1, ILevelRK4 - 1), TEMP3)
      Do I = 1, NTOT3T - NPWALL
        X1(I) = oX1(I)
     &         + dTempStep(ILevelRK4) * DT * TEMP1(I)
        If (X1(I).GT.XMAX) X1(I) = X1(I) - CUBESIZE
        If (X1(I).LT.XMIN) X1(I) = X1(I) + CUBESIZE
        If (DIMM.EQ.DIMM3D) Then
          X2(I) = oX2(I)
     &         + dTempStep(ILevelRK4) * DT * TEMP2(I)
          If (X2(I).GT.XMAX) X2(I) = X2(I) - CUBESIZE
          If (X2(I).LT.XMIN) X2(I) = X2(I) + CUBESIZE
        Endif
        X3(I) = oX3(I) 
     &         + dTempStep(ILevelRK4) * DT * TEMP3(I)
        If (X3(I).GT.XMAX) X3(I) = X3(I) - CUBESIZE
        If (X3(I).LT.XMIN) X3(I) = X3(I) + CUBESIZE
        If (tagRegroup .NE. 0) Then
          Call Regroup
        Endif
      Enddo
      Endif
      Endif  ! ALE
      Endif  ! ILevelRK4


      !If (Shepard .EQ. 1) Then
        Call Calc_ShepardCoef
      !Endif

      If (ifDual.EQ.1) Then
        Do I = 1, GridDual2
          Call Interpolation(X1Dual(I), 0.0, X3Dual(I),
     &                       Rho, RhoDual(I))
          Call Interpolation(X1Dual(I), 0.0, X3Dual(I),
     &                       U1, U1Dual(I))
          Call Interpolation(X1Dual(I), 0.0, X3Dual(I),
     &                       U3, U3Dual(I))
        Enddo
        Do I = 1, NTOT3T
          Call IntepGrad(X1(I), X3(I), U1Dual, DU1DX(I), DU1DZ(I))
          Call IntepGrad(X1(I), X3(I), U3Dual, DU3DX(I), DU3DZ(I))
        Enddo
      Endif

      If (ifDual.EQ.1) Then
        Call Calc_dRho_dT_Dual
      Else
        Call Calc_dRho_dT
      Endif

      Call PressureField

      Do I = 1, NTOT3T
        kRho(I, ILevelRK4) = Temp1(I)
      Enddo

      If (ifDual.EQ.1) Then
        Do I = 1, GridDual2
          Call Interpolation(X1Dual(I), 0.0, X3Dual(I),
     &                       P, PDual(I))
        Enddo
        Do I = 1, NTOT3T
          Call IntepGrad(X1(I), X3(I), PDual, DPDX(I), DPDZ(I))
        Enddo
      Endif

      If (ifDual.EQ.1) Then
        Call Calc_dU_dT_Dual
      Else
        Call Calc_dU_dT
      Endif

      Do I = 1, NTOT3T
        kU1(I, ILevelRK4) = Temp1(I)
        If (DIMM.EQ.DIMM3D) kU2(I, ILevelRK4) = Temp2(I)
        kU3(I, ILevelRK4) = Temp3(I) + Gravity
      Enddo

      If (ALE .EQ. ALEEuler) Then
        Call Grad(U1, DU1DX, DU1DY, DU1DZ)
        If (DIMM.EQ.DIMM3D) Call Grad(U2, DU2DX, DU2DY, DU2DZ)
        Call Grad(U3, DU3DX, DU3DY, DU3DZ)
        Call Calc_NLTerm
        Do I = 1, NTOT3T
          kU1(I, ILevelRK4) = kU1(I, ILevelRK4) + Temp1(I)
          If (DIMM.EQ.DIMM3D) kU2(I, ILevelRK4) =
     &                        kU2(I, ILevelRK4) + Temp2(I)
          kU3(I, ILevelRK4) = kU3(I, ILevelRK4) + Temp3(I)
        Enddo
      Endif

      If (NU .NE. 0.0) Then
      Call Calc_VisTerm
      Do I = 1, NTOT3T
        kU1(I, ILevelRK4) = kU1(I, ILevelRK4) + Temp1(I)
        If (DIMM.EQ.DIMM3D) kU2(I, ILevelRK4) =
     &                      kU2(I, ILevelRK4) + Temp2(I)
        kU3(I, ILevelRK4) = kU3(I, ILevelRK4) + Temp3(I)
      Enddo
      Endif

      If (ALE .EQ. ALELagrange) Then
        Do I = 1, NTOT3T
          kX1(I, ILevelRK4) = U1(I) 
     &       + dTempStep(ILevelRK4) * DT * kU1(I, ILevelRK4)
          kX2(I, ILevelRK4) = U2(I)
     &       + dTempStep(ILevelRK4) * DT * kU2(I, ILevelRK4)
          kX3(I, ILevelRK4) = U3(I)
     &       + dTempStep(ILevelRK4) * DT * kU3(I, ILevelRK4)
        Enddo
      Endif

      Enddo ! ILevelRK4


      Do I = 1, NTOT3T
        Rho(I) = oRho(I)
     &         + DT / 6.0 * kRho(I, 1)
     &         + DT / 3.0 * kRho(I, 2)
     &         + DT / 3.0 * kRho(I, 3)
     &         + DT / 6.0 * kRho(I, 4)
        U1(I) = oU1(I) 
     &         + DT / 6.0 * kU1(I, 1)
     &         + DT / 3.0 * kU1(I, 2)
     &         + DT / 3.0 * kU1(I, 3)
     &         + DT / 6.0 * kU1(I, 4)
        U2(I) = oU2(I) 
     &         + DT / 6.0 * kU2(I, 1)
     &         + DT / 3.0 * kU2(I, 2)
     &         + DT / 3.0 * kU2(I, 3)
     &         + DT / 6.0 * kU2(I, 4)
        U3(I) = oU3(I) 
     &         + DT / 6.0 * kU3(I, 1)
     &         + DT / 3.0 * kU3(I, 2)
     &         + DT / 3.0 * kU3(I, 3)
     &         + DT / 6.0 * kU3(I, 4)
      Enddo

      If (ALE .EQ. ALELagrange .AND. tagAVERMOVE .EQ. 1) Then
        Do I = 1, NTOT3T
          kX1(I, 4) = 
     &         1.0 / 6.0 * kX1(I, 1)
     &         + 1.0 / 3.0 * kX1(I, 2)
     &         + 1.0 / 3.0 * kX1(I, 3)
     &         + 1.0 / 6.0 * kX1(I, 4)
          kX2(I, 4) = 
     &         1.0 / 6.0 * kX2(I, 1)
     &         + 1.0 / 3.0 * kX2(I, 2)
     &         + 1.0 / 3.0 * kX2(I, 3)
     &         + 1.0 / 6.0 * kX2(I, 4)
          kX3(I, 4) = 
     &         1.0 / 6.0 * kX3(I, 1)
     &         + 1.0 / 3.0 * kX3(I, 2)
     &         + 1.0 / 3.0 * kX3(I, 3)
     &         + 1.0 / 6.0 * kX3(I, 4)
        Enddo
        CALL AVER(kX1(1, 4), TEMP1)
        CALL AVER(kX2(1, 4), TEMP2)
        CALL AVER(kX3(1, 4), TEMP3)
      Endif

      Do I = 1, NTOT3T - NPWALL
      If (ALE .EQ. ALELagrange) Then
      If (tagAVERMOVE .EQ. 0) Then
        X1(I) = oX1(I) 
     &         + DT / 6.0 * kX1(I, 1)
     &         + DT / 3.0 * kX1(I, 2)
     &         + DT / 3.0 * kX1(I, 3)
     &         + DT / 6.0 * kX1(I, 4)
        If (X1(I).GT.XMAX) X1(I) = X1(I) - CUBESIZE
        If (X1(I).LT.XMIN) X1(I) = X1(I) + CUBESIZE
        If (DIMM.EQ.DIMM3D) Then
          X2(I) = oX2(I) 
     &         + DT / 6.0 * kX2(I, 1)
     &         + DT / 3.0 * kX2(I, 2)
     &         + DT / 3.0 * kX2(I, 3)
     &         + DT / 6.0 * kX2(I, 4)
          If (X2(I).GT.XMAX) X2(I) = X2(I) - CUBESIZE
          If (X2(I).LT.XMIN) X2(I) = X2(I) + CUBESIZE
        Endif
        X3(I) = oX3(I) 
     &         + DT / 6.0 * kX3(I, 1)
     &         + DT / 3.0 * kX3(I, 2)
     &         + DT / 3.0 * kX3(I, 3)
     &         + DT / 6.0 * kX3(I, 4)
        If (X3(I).GT.XMAX) X3(I) = X3(I) - CUBESIZE
        If (X3(I).LT.XMIN) X3(I) = X3(I) + CUBESIZE
      Else ! tagAVERMOVE
        X1(I) = oX1(I) 
     &         + DT * TEMP1(I)
        If (X1(I).GT.XMAX) X1(I) = X1(I) - CUBESIZE
        If (X1(I).LT.XMIN) X1(I) = X1(I) + CUBESIZE
        If (DIMM.EQ.DIMM3D) Then
          X2(I) = oX2(I) 
     &         + DT * TEMP2(I)
          If (X2(I).GT.XMAX) X2(I) = X2(I) - CUBESIZE
          If (X2(I).LT.XMIN) X2(I) = X2(I) + CUBESIZE
        Endif
        X3(I) = oX3(I) 
     &         + DT * TEMP3(I)
        If (X3(I).GT.XMAX) X3(I) = X3(I) - CUBESIZE
        If (X3(I).LT.XMIN) X3(I) = X3(I) + CUBESIZE
      Endif ! tagAVERMOVE

      Endif
      Enddo

        If (tagRegroup .NE. 0) Then
          Call Regroup
        Endif

      Call PressureField

      End Subroutine
